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Q ; We report on the angular correlation function of Lyman-break galaxies 

(LBGs) at z ~ 4 and ~ 5 from deep samples obtained from the Great Ob- 
servatories Deep Origins Survey (GOODS). As for LBGs at z ~ 3, the shape of 
w(0) of the GOODS LBGs is well approximated by a power-law with the slope 
O ' f3 ~ 0.6 at angular separation 6 > 10 arcsec. The clustering strength of z ~ 4, 
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with brighter galaxies more strongly clustered than fainter ones. At smaller sepa- 
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rations, w(8) significantly exceeds the extrapolation of the large-scale power-law 



fit, implying enhanced spatial clustering on comoving scales r < 1 Mpc. We 
also find that bright LBGs statistically have more faint companions on scales 
6 < 20 arcsec than fainter ones. The enhanced small-scale clustering is very 
likely due to sub-structure, reflecting multiple galaxies within the same massive 
halos. A simple model for the halo occupation distribution along with the halo 
mass function in a ACDM cosmology, reproduces well the observed w(9). The 
scaling relationship of the clustering strength with volume density and with red- 
shift is quantitatively consistent with that of CDM halos. If we associate LBGs 
with dark matter halos that have the same clustering strength, this luminosity 
dependence of w(6) implies a close correlation between the halo mass and the 
star formation rate. A comparison of the clustering strength of three samples of 
equal luminosity limit at z ~ 3, 4 and 5 shows that the LBGs at z ~ 5 are hosted 
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in halos about 5-10 times less massive than those at z ~ 3 - 4, suggesting that 
star-formation was more efficient at z ~ 5. 

Subject headings: cosmology: observations — galaxies: distances and redshifts 
- galaxies: evolution — galaxies: formation 

1. Introduction 

In the current theoretical framework of galaxy evolution, the observed spatial clustering 
of galaxies is interpreted as reflecting the clustering of the dark matter halos that host 
them. Galaxies of different luminosity, spectral type, morphology, etc. are observed to have, 
in general, different clustering properties, because the specific criteria adopted to select 
them also select different "types" of halos. For example, because more massive halos have 
stronger spatial clustering than less massive ones, galaxies selected by criteria that also select 
the high-end of the halo mass function will be observed to be more clustered than galaxies 
whose selection criteria correspond to less massive halos. Similarly, the apparent "evolution" 
of the spatial clustering of galaxies as a function of redshift is interpreted as the result of 
the specific "mapping" of the selection criteria of the galaxies into a mass range of the halo 
mass function and of the evolution of the clustering of the halos, which is driven by gravity. 

This dependence of galaxy clustering on the selection criteria and how they map into 
properties of the halos is at the same time bad and good news. It is bad news because it 
makes it very difficult to use galaxies to infer anything about the evolution of dark matter, 
and hence about cosmology, unless one knows exactly how to go from galaxies' properties to 
halos' properties, including taking into account the "astrophysical" evolution of the former, 
a very complex problem. It is good news because, at least for certain type of galaxies and in 
some simple cases, one can expect to learn something abut the physical relationship between 
certain properties of the galaxies and those of the underlying halos, possibly testing some of 
the key ideas behind our understanding of galaxy evolution. 

One important case where the relationship between galaxies and halos is being success- 
fully explored by means of the analysis of the spatial clustering is that of the Lyman-break 
galaxies. These are star-forming galaxies whose rest-frame UV spectral energy distribution 
is not severely obscured by dust. This type of source is very effectively selected by means 
of their broad-band UV colors, constructed from around the 912 A Lyman limit, the Ly« 
forest region (pronounced in broad-band photometry at z > 2), and the otherwise relatively 
featureless continuum redward of the Lya line (e.g. Steidel & Hamilton 1993; Steidel et al. 
1995; Madau et al. 1995). Using optical passbands, where one can take advantage of the 
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low sky and of the sensitivity and large area coverage of CCD imagers and spectrographs, 
the technique is effective at redshifts 2.5 < z < 6. Large surveys of Lyman-break galaxies 
have been conducted over the past ten years in this redshift range, resulting in large and 
well controlled samples, including large spectroscopic ones (e.g. Steidel et al. 1999; 2003; 
Giavalisco et al. 2004; Dickinson et al. 2004; Sawicki & Thompson 2005; Vanzella et al. 
2005; E. Vanzella et al., in preparation). These surveys are also well suited for clustering 
measures. 

Lyman-break galaxies have been found to have relatively strong spatial clustering, with 
a correlation length comparable to that of present-day bright spiral galaxies (Giavalisco et 
al. 1998; Giavalisco & Dickinson 2001; Porciani & Giavalisco 2002; Adelberger et al. 2005). 
This has been interpreted as evidence that they are biased tracers of the mass density field 
(Giavalisco et al. 1998; Adelberger et al. 1998; Giavalisco & Dickinson 2001; Bagla 1998). 
This means that these galaxies are preferentially hosted in the most massive dark-matter 
halos (Giavalisco & Dickinson 2001; Foucaud et al. 2003; Adelberger et al. 2005), which 
have stronger spatial clustering than less massive halos (see, e.g. Mo & White 1996). A key 
piece of evidence linking the mass of the halo to physical properties of the galaxies came 
with the discovery that the strength of the spatial clustering of Lyman-break galaxies is 
a function of their UV luminosity, with brighter galaxies having larger correlation length 
(Giavalisco & Dickinson 2001; Foucaud et al. 2003; Ouchi et al. 2004; Adelberger et al. 
2005; Allen et al. 2005). In the halo interpretation, this is explained if brighter galaxies are, 
on average, hosted in more massive halos. Since the halo mass function is relatively steep, 
this also implies that the statistical relationship between mass and UV luminosity must be 
characterized by a correspondingly small scatter (Giavalisco & Dickinson 2001). Because the 
UV luminosity of Lyman-break galaxies is powered by star formation, this in turn means 
that the mass and the star-formation rate are relatively tightly correlated or, in other words, 
that the halo mass is a primary parameter that controls the star-formation activity of the 
Lyman-break galaxies as expected in galaxy formation models. 

The next logical step is to try to derive the relationship between mass and star-formation 
rate, thus testing the basic ideas behind the physical models of star-formation in young 
galaxies (e.g. White & Rees 1978). Furthermore, the dependence of the clustering strength 
on the UV luminosity (clustering segregation hereafter) contains direct information on the 
scaling relationship between the halo clustering strength (bias) and the volume density once 
the galaxy luminosity function is known. Since the shape of this relationship depends on 
the shape of the power spectrum, the observed clustering segregation can be used to test 
key predictions of the CDM theory (Giavalisco & Dickinson 2001). To do this, however, 
information on the halo sub-structure, or the halo occupation distribution function (HOD), 
is also required (e.g. Bullock et al. 2001, Berlind & Weinberg 2002). This function describes 
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the likelihood that one halo contains more than one galaxy and how this depends on its mass. 
Once the HOD and the clustering segregation have been measured, then the relationship 
between mass and UV luminosity (star-formation rate) can be constrained, and a successful 
model will have to simultaneously reproduce the observed clustering segregation and UV 
luminosity function. 

In principle, the presence of sub-structure and the HOD function can be constrained 
by the shape of the angular correlation function and, as we will see later, by the statistics of 
close pairs of galaxies. In practice, this requires measures with relatively high S/N from large 
and faint samples, unavailable until recently, since most of the secondary halo-companion 
galaxies have lower mass than the central galaxy and are thus less luminous on average. 
Furthermore, large samples covering a sufficiently wide dynamic range in luminosity are 
necessary to provide a robust measure of the shape of the clustering segregation. Giavalisco 
& Dickinson (2001) first reported the detection of the clustering segregation and measured 
its shape but with low S/N. Using a deep sample of LBGs at z ~ 4 and 5 obtained with 
the Subaru telescope and the Suprime camera, Hamana et al. (2004) reported evidence 
that the angular correlation function of LBGs deviates from a power-law at small angular 
scales and derived the HOD function to explain the deviation. Using the largest sample 
of LBGs to date at z ~ 3, Adelberger et al. (2005) measured the shape of the angular 
correlation function at large angular separations with high S/N. Their sample, however, is 
not deep enough to include halo companion galaxies in appreciable number and detect the 
effect of sub-structures; they could only set an upper limit of 5% to the number of galaxies 
in common halos. 

In this paper we present a study of the clustering properties of Lyman-break galaxies at 
z ~ 4 and ~ 5 at faint flux levels. We use two deep and relatively large samples obtained with 
EST and ACS during the program of multi-band imaging of the Great Observatories Origins 
Deep Survey (GOODS). In fact, these ACS samples of B 435 and V606~band dropouts are 
currently the largest and the most complete at their depth of such galaxies 1 . Since their size 
is suitable to measure the angular clustering with high S/N, we have embarked on a project 
to characterize the physical association between the activity of star formation of the galaxies 
and the properties of the dark matter halos, primarily the mass, thus providing quantitative 
tests to the CDM hierarchical theory of galaxy formation. Here we present significantly 
improved measures of the clustering segregation down to faint flux limits (-2850~ 27), show 



1 Analogous samples have also been obtained during the Hubble Ultra Deep Field survey, which covers 
a small portion of the GOODS southern field (CDF-S) with identical passbands as the primary GOODS 
survey. While significantly deeper than the GOODS samples, the HUDF ones are much smaller in size 
because of the minimal areal coverage of the HUDF. 
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direct evidence of sub-structure in halos, and derive the corresponding HOD. We postpone 
to a following paper the derivation of the relationship between mass and star-formation rate. 
All magnitudes in this paper are in the AB scale of Gunn & Oke (1975) and, when necessary, 
we use a cosmology with Q m = 0.3, Q\ = 0.7, a 8 = 0.9, T = 0.21, H = 100/i km s _1 Mpc -1 
with h = 0.7 and baryonic density Qb = 0.04. 

2. The Data and The Samples 

Deep optical multi-band imaging with HST and ACS has been acquired as part of the 
program of observations of the GOODS. These ACS data consist of a mosaic of images 
in each of the -B435, VJ306, ^775 and z 850 passbands for a total of 3, 2.5, 2.5 and 5 orbits of 
integration time, respectively, in each of the two GOODS fields. These two fields are centered 
around the Hubble Deep Field North (HDF-N) and Chandra Deep Field South (CDF-S) 
and cover approximately 10'xl7' each, for a total area coverage of about 0.1 square degree. 
To search for supernovae at high redshift (Riess et al. 2004), the observations in the V 60 6, 
2775 and ^850 bands were not done continuously, but sequenced into "epochs" separated by 
about 45 days from each other. Because the nominal roll of HST changes at a rate of about 
1 degree per day, this observing strategy resulted in images acquired during different epochs 
having position angle in the sky differing by integer multiple of 45 degree. This gives the 
characteristic jagged edges of the GOODS mosaics in the Vqqq, ^775 and ^§50 bands. Giavalisco 
et al. (2004a) present an overall description of the GOODS project, and we refer the reader 
to that paper for additional details on the ACS observations. 

The Vqqq, "^775 an d £850~band mosaics presented in Giavalisco et al. (2004a,b) only 
include the first 3 epochs of observations stacked together, while the samples of Lyman- 
break galaxies at z ~ 4 and ~ 5 that we are about to discuss here have been extracted from 
the version vl.O reduction of the ACS data, which include the full stack of the five epochs. 
The 5 435 -band mosaic discussed by Giavalisco et al. (2004a), however, already includes all 
the data, since all the £>-band observations were taken continuously. A complete description 
of the vl.O ACS images will be presented in an upcoming paper (M. Giavalisco et al., in 
preparation). Detailed information can also be found in the GOODS website, at the link 
http : //www . stsci . edu/science/goods/vl_release_readme/h_goods_vl . CLrdm . html. 

We now provide a brief summary of the data reduction and stacking process for the 
version vl.O data release. Differently from the version v0.5 public release and from the mosaic 
discussed in Giavalisco et al. (2004a), for the version vl.O reduction we have re-processed 
all the data (including the 5435-band images) using the best calibration and reference files 
available. We have improved treatments of the geometrical distortions of ACS using updated 
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coefficients together with corrections for the velocity aberration distortion, and we have also 
improved rejection of cosmic rays and other blemishes. We have stacked the individual 
exposures together in each band using the drizzling algorithm in two independent phases. 
In the first phase, images taken in the same filters are sky-subtracted and drizzled onto 
a common pixel grid at the instrument native scale (0.05 arcsec/pixel). Cosmic rays and 
deviant pixels are identified during this process and flagged in mask files specifically created 
for this purpose. Information on detector blemishes (hot pixels, bad columns, etc.) in post- 
pipeline masks (which have also been drizzled onto the same grid) is included in the new 
masks at this time. During the second phase, the images and the mask files are blotted 
back to the original positions, drizzled again onto a common astrometric grid with scale 0.03 
arcsec/pixel, and stacked together. During this process corrections for the ACS geometrical 
distortion are applied, cosmic rays flagged during the previous processing block are masked 
out from the stack, and additional, low-level cosmic rays and defects are identified and 
masked, too. The final stacks reach similar 1-er surface brightness limits for the north and 
south, which are listed in Table 1. 

Multi-band source catalogs, which are public and available at the Web site listed above, 
have been created using the SExtractor package (Bertin & Arnouts 1996). We made the 
detections using the z 850 -band mosaic, and then used a variety of photometric apertures, 
including isophotal and a suite of circular apertures of varying radius, as "fixed apertures" 
to carry out matched photometry in all bands. We have measured the completeness of the 
catalogs as a function of the size and magnitude of the sources through extensive Monte Carlo 
simulations. In these simulations we have inserted a very large number of model galaxies 
(169,100) into the ACS mosaic, after convolution with the ACS PSF and after adding the 
appropriate Poisson and detector noise. We used galaxies with both exponential and r 1 / 4 light 
profiles with random ellipticity and orientation in the sky, with varying apparent magnitude 
and size, and then retrieved them with identical procedures as for the real galaxies. From 
these simulations we estimate that the catalogs are ~ 50% complete at z 850 ~ 26.5 for sources 
with half-light radius ri/ 2 < 0.2 arcsec, and about 10% complete at £850~ 27. We have also 
run multi-band simulations for sources with a given input spectral energy distribution (SED) 
to test completeness relative to color selection, as we will discuss in a moment. 

We have also used deep ground-based multi-wavelength imaging data to select samples 
of Lyman-break galaxies at redshifts z ~ 3 ([/-band dropouts) in the CDF-S field, albeit 
in a larger area than the ACS fields. These data were acquired as part of the program of 
ancillary observations of the GOODS project. They are also described in Giavalisco et al. 
(2004a), and we refer the reader to that publication for details. Here we will summarize the 
key properties of these data. The B and R bands were taken with the Wide Field Imager 
(WFI) at the 2.2m MPG/ESO telescope, while the [/-band observations were carried out 
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with the MOSAIC II Camera at the CTIO 4m telescope. The stacked image in each band 
covers a contiguous region of approximately 0.4 degree 2 . Some parameters of these ground- 
based data sets are detailed in Table 2. 

Samples of star-forming galaxies at z ~ 3, 4 and 5 were photometrically selected using 
the Lyman Break technique (e.g. Steidel et al. 1999; Madau et al. 1996; Giavalisco et al. 
2004b.; see also Giavalisco 2002 for a review). The technique and its application to clustering 
studies (Giavalisco et al. 1998; Adelberger et al. 1998; Giavalisco & Dickinson 2001) have 
been extensively discussed in the past, and we refer the reader to the cited literature. 

Lyman-break galaxies (LBGs) at z ~ 4 and ~ 5, -B435 and V606~band dropouts, respec- 
tively, were extracted from the GOODS vl.O rl.l catalogs using the same color selection 
criteria described by Giavalisco et al. (2004b), which we report here for convenience. We 
defined the -B435 and V(j06 _ band dropouts as 

(.8435 - V m6 ) > 1.2 + 1.4 x (\/ 606 - ^850) A 

(-B435 — V e0e ) > 1-2 A 1^606 — ^850 < 1-2, 

and 



(^606 


- ^775) 
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1.5 + 0.9 x (i 775 -z 850 ) V 


(^606 


- ^775) 
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2.0 A (V(306 _ ^775) > 1-2 A 


(^775 


— ^850 ) 


< 


1.3, 



where the symbols V and A are the logical "OR" and "AND" operators, respectively. In both 
cases we have limited the apparent magnitude of the samples to z$5o< 27.0 for completeness. 
We have also visually inspected each galaxy included in the samples and removed objects 
with stellar morphology and obvious spurious detections, such as diffraction spikes. Down 
to the adopted flux limit of -2 850 < 27.0, the ACS samples include 2463 5 435 -band dropouts 
and 878 V606~band dropouts. 

Lyman-break galaxies at z ~ 3 were selected using 

(U-B) > (B — R) + 0.6 A 
{U-B) > 0.9 A (B - R) < 2.5 

This ground-based sample includes 1609 galaxies in the central region of the field which was 
used for the clustering analysis down to a flux limit of R < 25.5. The corresponding surface 
density values of our samples down to the above mentioned flux limits are £ = 2.7, 7.7, 2.7 
galaxies/arcmin 2 for the U, B i35 and V^-dropouts, respectively. 
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Contamination from galactic stars is nearly negligible for the -B435 and V606~band dropouts 
because they were morphologically culled taking advantage of the high angular resolution of 
the GOODS ACS data. For the ground-based sample such culling is not possible; however 
stellar contamination in high galactic latitude fields such as GOODS has been found to be 
~ 4% in previous spectroscopic surveys (Steidel et al. 1998). We also impose a minimum 
signal to noise ratio (S/N > 10) in order to filter out spurious sources. This also helps 
(relatively) homogeneous detection of star-forming galaxies up to a certain magnitude limit 
when dealing with data produced with several separate pointings with overlaps. Even at the 
cost of losing some fraction of the faintest LBGs, it is desirable for clustering studies to avoid 
an artificial clustering signature induced by interlopers and inhomogeneity of the data. 



3. Simulations 

At z ~ 3 the selection of Lyman-break galaxies ([/-band dropouts) has been charac- 
terized in great detail thanks to the systematic spectroscopic identification of thousands of 
candidates (e.g. Steidel et al. 2003). In particular, this extensive body of work has shown 
that the Lyman-break technique is very efficient, with a relatively low rate of contamination 
from low-redshift interlopers, and has yielded the redshift distribution function N(z) asso- 
ciated with a given set of color selection criteria. While spectroscopic identifications have 
been made of Lyman-break galaxies at z ~ 4 and ~ 5 (Steidel et al. 1999; Giavalisco et 
al. 2004), including a number presented in this work (see later), the spectroscopic samples 
still remain too small for accurate measures of the redshift distribution function associated 
with a specific color selection. Fortunately, as the z ~ 3 case has demonstrated (Steidel et 
al. 1999), Monte Carlo simulations are very effective in providing robust estimates of N(z), 
which are needed to derive the real-space clustering from the angular one and to measure 
the galaxies' volume density. 

The simulations consists of creating artificial LBGs over a large redshift range (we used 
2 < z < 8) with assumed distribution functions of UV luminosity, SED, morphology and 
size, inserting them in the real images with random orientations and inclinations and with 
the appropriate amount of photon noise, detecting them, and measuring their "observed" 
photometry and morphology as if they were real galaxies. We have adjusted the input 
distribution functions of UV luminosity (UV luminosity function), size and SED so that the 
observed distribution functions of apparent magnitude, UV colors and size of the simulated 
galaxies best reproduce the same quantities of the real galaxies. Specifically, we used galaxies 
with exponential and r 1 / 4 light profile in equal proportions and size extracted from a log- 
normal distribution function, as described by Ferguson et al. (2004). This method was first 
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used by Giavalisco et al. (2004b) and Ferguson et al. (2004); its application to the measure 
of the UV luminosity function of LBGs will be further discussed in an upcoming paper (M. 
Giavalisco et al., in preparation). 

The average redshifts obtained from the simulations are z— 3.2 ± 0.3, 3.8 ± 0.3, 4.9 ± 0.3 
for the U, -B435 and V606~band dropouts, respectively. Figure 1 shows the redshift distribution 
N(z) of simulated LBGs together with that of a subset of LBGs with spectroscopic redshifts. 

4. Measuring the clustering properties of Lyman— break Galaxies 

We now discuss the measure of angular and real-space clustering of LBGs and other 
clustering properties that we have carried out from our samples. The primary measure is 
that of the angular correlation function w{6), which we deproject with the Limber transform 
using the redshift distribution function estimated from the simulations (but see later for a 
discussion on our spectroscopic observations) to derive the real-space correlation length, in 
the power-law approximation. In our faint samples, the angular correlation function is very 
well approximated by a power law at large angular separations, roughly 6 > 10 arcsec, but 
it exceeds the extrapolation of the large-scale power-law fit at small scales. To understand 
the nature of this excess clustering, which is not observed in bright samples, either our own 
[/-band dropout sample or even larger ones (e.g. Adelberger et al. 2005), we have carried 
out an analysis of the mean number of close neighbors around both bright and faint galaxies. 
We have used the halo occupation distribution (HOD) formalism to predict the shape of w(9) 
and compare it with our measure. 

4.1. The angular correlation function 

The inversion of the angular correlation function is a robust method to derive the real- 
space correlation length if the redshift distribution function is well known (see Giavalisco 
et al. 1998; Giavalisco & Dickinson 2001; Adelberger et al. 2005). The measure of w(9) 
is relatively straightforward; however some care is required in the analysis of the random 
errors, and, above all, of systematic errors, since they can significantly bias the result. Given 
the relatively small area covered by the GOODS fields and the way the observations have 
been obtained, the two most significant sources of systematic errors are cosmic variance (the 
integral constraint bias) and sensitivity fluctuation in the survey, which can mimic a spurious 
clustering signal. These will be discussed separately in dedicated subsections. 

We have estimated the angular correlation function using the estimator proposed by 



-10- 



Landy and Szalay (1993): 



DD(6)-2DR(6) + RR(6) 

Wobs{6) = rW) ' (1) 



where DD(6) is the number of pairs of observed galaxies with angular separations in the 
range (6 — 56/2,6 + 89/2), RR(6) is the analogous quantity for random catalogs with the 
same geometry as the observed catalog and DR{6) is the number of random galaxy cross 
pairs. We have chosen logarithmic binning of the angular separations to provide adequate 
sampling at small scales while avoiding excessively fine sampling at large scales which would 
have resulted had linear binning been adopted. We have repeated the measures using 4 
different (logarithmic) binsizes for each sample. 

To account for the random errors, we have estimated the error bars for each angular bin 
by bootstrap resampling of the data (e.g. Benoist et al. 1996, Ling, Barrow & Frenk 1986). 
We have computed values of w(6) from 100 randomly selected subsets of galaxies and then 
used the standard deviation of w{6) as our primary error estimate. This method, however, 
only accounts for random errors appropriate to our sample size, and does not include the 
uncertainty due to the finite size of our samples, because we use the observed galaxies for 
the error estimate. In other words, the bootstrap method fails to account for fluctuations 
on scales larger than our survey. We have used Monte Carlo simulations to estimate these 
fluctuations of w{6) due to cosmic variance and validate our accounting of random errors. 
We have created a large number of realizations (50) of GOODS samples at z ~ 4 and ~ 5 
extracted from a population with the same large-scale clustering properties (i.e. surface 
density and angular correlation function — see Porciani & Giavalisco 2002 for a description 
of the method) as the observed galaxies. Since the goal of this test is to estimate the 
uncertainty induced in the observed w{6) by the large-scale fluctuations of the density field 
of the parent population of galaxies from which are samples are extracted from, neglecting 
the small-scale clustering is inconsequential. Each realization simulates an area of the sky 
of 1 degree 2 within which we have selected the mock "GOODS" samples of LBGs. We have 
measured w(6) in each sample, and derived the standard deviation of the 50 measures at 
each angular bin and used it as an estimate of the error bar, finding consistent results with 
the bootstrap method. 

To test the dependence of the angular clustering on the luminosity of the galaxies, we 
have measured the correlation function for various magnitude cuts in each sample. For the 
-B435 and V606~band dropout samples we have defined sub-samples at z$5o< 27.0, 26.5 and 
26.0; for the [/-band dropouts we have used R < 25.5, 24.5 and 24.0. In what follows, we 
will denote each sub-sample as B(V)270, 265, 260 and U255, 245, 240, respectively. The 
number of galaxies included in each sub-sample is 2463, 1410 and 858 for the -B 435 -band 



dropouts, 878, 407 and 218 for the V^-band dropouts, and 1609, 1050 and 554 for the 
[/-band dropouts. 



4.2. Correction for sensitivity variations 

A possible source of systematic error in surveys such as GOODS, which consist of either 
mosaics of images (the ACS images) or images taken with CCD mosaics (the ground-based 
images), is fluctuations of sensitivity across the surveyed area. Such fluctuations, which 
typically occur on scales that are a fraction of the linear size of the survey, can be sufficiently 
strong to mimic the presence of clustering and bias the measures of w{6). 

To estimate this effect in LBGs samples is not straightforward, because these are color- 
selected as well as flux-limited, and thus the final selection function is a nontrivial function 
of sensitivity variations in all the different passbands used in the color selection as well 
as of their respective limiting magnitudes. The best way to take all of these factors into 
account is to use Monte Carlo simulations. We have generated over 100,000 artificial U 
and 5 435 -band dropouts and placed them at random positions in each of the U, B, R and 
-B435, ^606, 2775, -^850 images of the ground-based and ACS survey, respectively (see above for 
a description of the simulation technique). We have subsequently "detected" the artificial 
galaxies and subjected them to the same color selection as if they were real U or 5 435 -band 
dropouts. To the extent that this process is a faithful representation of the overall selection 
function of the survey, the spatial distribution of these color-selected mock galaxies will 
reflect various systematics including the sensitivity variations across the field, the effects of 
holes due to saturated sources or large seeing (of our ground-based mosaics), to the applied 
color-selection and so forth. 

The correction factor can be estimated by measuring the angular correlation function, 
w s im{0), against a purely random distribution with the same geometry (this should be essen- 
tially zero for perfectly homogeneous data). Since w(8) is the excess probability of finding a 
pair, our previous measurement can be corrected as 

1 . ( a\ l+W raw (9) 

1 + w corr {e) = — — (2) 

l + w sim (9) 

where w raw is the measurement using the LS estimator, w sim is the correlation function of 
artificial galaxies and w corr is the corrected angular correlation function of the dropouts. 
While the correction is non-negligible at small scales for the [/-dropouts, it has negligible 
effect on the ACS data, hence the correction is applied only to the £7-dropout measurement. 
The subscript of w corr {6) will be dropped hereafter unless stated otherwise. 
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4.3. Integral Constraint 

Due to the finite size of any given survey, the correlation function measured from a 
sample is underestimated by a constant known as the integral constraint, referred to as IC 
hereafter: 

(9) = u) obs (6) + IC (3) 

Since the measure of w{6) is underestimated by the same amount on all scales, this 
also leads to an overestimate of the correlation slope, (3. This can be estimated by doubly 
integrating u) true over the survey area (Roche & Earles 1999). 

ic = ^ J Ju true (e)dn 1 n 2 

T li RR(9 i ) H i RR(9 i ) 

Note that the integral constraint depends on both the size of a given survey (Q) and 
the galaxy power spectrum (approximated by a power-law ; A w and (3). However, once 
the correlation slope (3 is fixed, IC/A W only depends on the size and the shape of the 
survey area. Therefore, one needs to make robust estimate of (3 from the observational 
data, then for a fixed f3, IC/A W can be estimated. We calculate the ICs from random 
catalogs generated over a field with the same size and the geometry as our survey area with 
a range of and (3. For each set (A w , (3) we calculate IC value using Equation (4) and 
compute the \ 2 °f the measured w f, s (8). Once f3 and IC/A W are robustly measured from 
the full sample, for other flux-limited samples, A w is the only parameter that needs to be 
fitted as w mo dei = A W (9~P — (IC/A w )q) provided that (3 does not change significantly with 
samples of different flux limits. This helps to make reliable measures of the amplitude A w 
for sub-samples that have lower signal-to-noise than the full sample. 

For the B270 sample, we find that a set (A w , (3, IC) = (0.31, 0.6, 0.012) best describes 
the measures, however, it should be noted that there is degeneracy between A w and f3. In 
other words, there is a range of acceptable A w s and /3s for a fixed IC value even for high S/N 
data, and hence one needs to somewhat arbitrarily fix the slope within the given range. We 
fix (3 = 0.6 and IC/A W = 0.039 for all B 435 and V^-band dropout samples. Similarly, for the 
[/-band dropouts, we find (A w , (3, IC) = (0.50, 0.6, 0.012) thus IC/A W = 0.024. Note that 
IC/A W is much smaller for the [/-band dropouts due to the larger field size. On the other 
hand, the IC values for the full sample of the -B435 and t/-band dropouts are comparable, 
because the faint ACS -B435 and V6O6 dropout samples are less clustered (i.e. smaller A w ) 
than the brighter ground-based sample. Using this method, we find IC = 0.012, 0.024, 0.031 
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for the U255, U245 and U240 sample, IC = 0.012, 0.020, 0.027 for the B270, B265 and B260 
sample and IC = 0.024, 0.031, 0.043 for the V270, V265 and V260 sample, respectively. 

When the slope is allowed to vary, the procedure becomes iterative because the derived 
IC depends on w(9) and vice versa, and one cannot conveniently set IC/A W to be a constant 
for the fitting. As Adelberger et al. (2005) pointed out, this is often very unstable regardless 
of the initial guess. Therefore, as a consistency check we estimated the IC using an alternate 
method suggested by Adelberger et al. (2005). Their method takes advantage of the fact that 
matter fluctuations, (Tq DM within a given survey volume are in the linear regime and can be 
estimated from the power spectrum with an appropriate window function that accounts for 
the geometry of the survey volume (see their Equation 20). By definition, IC, the variance 
of density fluctuation of galaxies is IC as b 2 o 2 CDM . We estimate the galaxy bias, b, from 
the correlation function itself by using the relation b = a 8t g/as(z), where ag, g is the galaxy 
variance in spheres of comoving radius 8 hr l Mpc (Peebles 1980, Equation [59.3]): 

2 72(r /8h- 1 Mpc)' 1 
a& ' 9 ~ (3- 7 )(4- 7 )(6- 7 )2T 

where r is the correlation length inferred from the angular correlation function (see next sub- 
section), 7 = /3 + 1 and a 8 (z) is the linear matter fluctuation on the same scale extrapolated 
from cr 8 (0) = 0.9. 

When using this method, we find values of both IC and (5 that are consistent with 
the previous method, but the IC values are slightly smaller in most cases. For example, 
for the B270 sample, we find IC = 0.010 instead of 0.012, for B265 0.018 instead of 0.020, 
which makes little difference in the fitted A w or (3 values. This may imply that the generic 
correlation function slope j3 may be slightly steeper than the fiducial value 0.6 for some 
sub-samples as can be seen in Table 3. For a fixed r , cr$ tg declines with the slope 7 at 
1.5 < 7 < 2.0. However, for a fixed 7 (or (3), both methods are in good agreement. 

Since the IC is u>(#)-dependent, power-law fitting and correcting for the IC are in fact 
simultaneous processes. In each iteration, only the w b s {9) measurement with 6 > 10 arcsec 
is used for the power-law fit (to probe large scale structure) because at 9 < 10 arcsec there 
may be significant contribution from highly nonlinear small-scale clustering. Figure 2 shows 
our results for the full -B 435 -band (left panel) and V6 6-band (right panel) dropout sample 
together with the best power-law fit. Clear departure from a pure power-law is apparent in 
both cases on small scales {6 < 10 arcsec). 

We fit the data to a linear function in logarithmic space. To estimate la confidence 
intervals for the parameters A w and (3, we carried out a large ensemble of random realizations 
of the measured w{6) assuming normal errors, and calculated best-fit parameter values for 



(5) 
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each of these synthetic data sets. Because we select randomly among the w(6) measures of 
different binsizes for each realization, the peak values of A w and j3 may be slightly different 
from those derived from the individual w(9) measures mentioned in the previous subsection. 
The fitted values of A w and f3 are provided in Table 3. 



4.4. The Real-Space Correlation Function £(r) 

We derive the real-space correlation function £(r) by inverting w{9) using the Limber 
transform (Peebles 1980). If the real-space correlation function has the form £(r) = (r/ro)~ 7 , 
the angular correlation function also has to be a power-law, w{9) = A w 9~@ where /3 = 7 — 1 
and A w is related to £(r) as, 

A w = Crl J F(z)D l -\z)N(z) 2 g(z)dz x [ J N(z)dz]- 2 

where D e is the angular diameter distance and N(z) is the redshift selection function derived 
from the simulations. 

g(z) = ^l(i + z) 2 (i + n z + n A ((i + z)- 2 -i)y/ 2 } 
r C ^ r[( 7 -i)/2] 
° = ^ r (7 / 2 ) • 

We have used the same N(z) for all the different limiting magnitudes since the simula- 
tions indicate little change in the redshift distribution for each case. Note that A w and (3 are 
obtained from the measures at 9 > 10 arcsec to be representative of large-scale clustering. 
It is interesting, however, that if we include the data points at 9 < 10 arcsec for the fit and 
allow (3 to vary, both A w and f3 get overestimated by an amount that changes tq very little. 
For example, for the B270 sample, we obtained the ro-value of 2.82 ±0.20 h~ l Mpc (2.81^q \ x 
h^ 1 Mpc) when using the data points at 9 > 3 (9 > 10) arcsec. Hence our estimates of ro 
are robust regardless of the fiducial value chosen to be considered as large-scale. 



The correlation lengths ro derived using this method are listed in Table 3. Two different 
tq values derived from the fits with and without the correlation slope j3 fixed to a fiducial 
value 0.60, are consistent with each other within the errors. In the following discussion, we 
regard the former as our best r values. The best-fit r values range from 2.8 h' 1 Mpc to 
7 h~ x Mpc (comoving) depending on the median luminosity of the sub-samples. In each 
redshift range, brighter samples have larger correlation lengths. We have compared our 
measures with those from other groups for samples with similar median luminosity and in 
the same redshift range generally finding good agreement. At z ~ 3 Adelberger et al. (2005) 
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found r = 4.0 ± 0.6 h 1 Mpc for galaxies with magnitude 23. 5< 1Z <25.5 (z = 2.9), in 
good agreement with our R < 25.5 sample for which we found r = 4.0 ± 0.2 h" 1 Mpc. At z 
~ 4 and 5, Ouchi et al. (2004) reported r = 4.1 ±0.2 and 5.9+{;^ h' 1 Mpc for their LBGMs 
(i f < 26) and LBG,25s (z' < 25.8) samples, respectively. Using an LBG template spectrum 
that reproduces the median UV colors of the samples, we computed that the flux limit of 
Ouchi et al. (2004) correspond to z 850 < 26.0 and 25.8, for which we find r = 5.1^5 and 
5. 30^;[;q hr l Mpc. While the correlation lengths of the two samples of V^-band dropouts 
are in good agreement, we noticed that Ouchi et al. reported a lower value than ours for the 
-E?435-band dropout sample, although the two measures overlap at the ~ 1.3a level. 



5. Results 

5.1. Dependence of Clustering strength on Luminosity 

A number of authors (Giavalisco & Dickinson 2002; Foucaud et al. 2003; Adelberger et 
al. 2005; Ouchi et al. 2004) have reported the existence of the clustering segregation with the 
UV luminosity for the LBGs, namely of the fact that the galaxies with higher UV luminosity 
have stronger spatial clustering (e.g. larger spatial correlation length) than fainter ones. We 
have examined the dependence of the clustering strength of our samples on luminosity using 
several apparent magnitude cuts constructed from the full samples. Apparent magnitude cuts 
roughly correspond to absolute magnitude cuts because the luminosity-distance changes only 
by rs 20% over the redshift range covered by each of our samples and the redshift selection 
functions are relatively narrow. Whether or not we fix the correlation slope to a fiducial 
value of P — 0.6 (Porciani & Giavalisco 2002), we find that the correlation amplitude A w , 
and ro accordingly, increases with median luminosity. 

The correlation length of the 5 435 -band dropouts increases from 2.8±0.2 hr l to 3.7±0.3 
hr x to 5.1^5 hr l Mpc with increasing luminosity of the samples (see Tables 3). 

Note that the faintest sample of _B 435 -band dropouts, which is our faintest sample at 
any redshift, reaching absolute magnitude M 1700 = —18.52, also has the smallest correlation 
length, r = 2.8 h~ l Mpc. A similar trend is also found for the V^-band dropouts. For 
a fixed absolute luminosity, objects will be fainter at z ~ 5 than their counterpart at z ~ 
4 by 0.6 - 0.7 mag in 2 85 o-band due to cosmological dimming. This suggests that the 
correlation lengths for the -B435 and V606~band dropouts for a fixed absolute UV luminosity 
are comparable (see Table 3 for details). 

In agreement with previous results, we find that the ?7-band dropouts also follow a 
similar trend. For the U255 sample (R < 25.5) we measure a correlation length r = 4.0lg;2 
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h 1 Mpc, in an excellent agreement with the value reported by Adelberger et al. (2005), 
4.0 h- 1 Mpc for roughly the same magnitude threshold ((R — TV) ~ 0.03 at z ~ 3 using 
our mean LBG spectrum). The correlation length increases to r = 7.8^o'g hr x for R < 
24.0. Interestingly, this is comparable to 8 h^ 1 Mpc, the correlation length of a sample of 
distant-red galaxies (DRGs) found at a slightly lower redshift range (Daddi et al. 2003). 

The left three panels of Figure 3 show w(8) for the full sample (solid lines) and the 
brightest sub-samples (dashed lines) for each dropout flavor together with best-fit power- 
laws. The power-fits when (3 is fixed to 0.6 are shown for clarity. The right panels show the 
derived correlation lengths for these samples. All samples show a clear trend that brighter 
samples are more strongly clustered in real space. 

The sub-samples that we have used to measure the clustering segregation are not in- 
dependent from one another, because the galaxies that comprise the brighter samples also 
belong to the faintest one thus diluting the weak clustering signal from the faintest galaxies. 
We have repeated the measures using two mutually independent samples that we built by 
splitting the full fi 435 -band dropouts sample into two independent magnitude bins, 26.3 < 
z 85 o< 27.0 and z^ Q < 26.3 with roughly equal number of galaxies in each bin. For these two 
samples we have measured r = 4.0 ± 0.3 and 2. 3^0^ h^ 1 Mpc, respectively. Note that the 
correlation length of the faint sample is now smaller than 2.8 ±0.2 Mpc of the z$5o< 27.0 
sample, consistent with the luminosity segregation. 

Recently, Norberg et al. (2001) and Zehavi et al. (2002) observed luminosity segrega- 
tion in the optically selected galaxies in local universe with high S/N. Note that a direct 
comparison between our results and the local ones is not possible, since the UV luminosity 
is powered by star formation rate, while the optical luminosity is mostly powered by stellar 
mass. However, understanding how UV luminosity scales with the clustering strength, at 
both high and low luminosities, will shed light on the relationship between halo mass and 
star formation rate. Unfortunately, our brightest magnitude cut made to the ACS sample 
is -2850 = 26.0, which is still more than 1 mag fainter than m* (m* is £850~ 24.7 and 24.8 at 
z = 4 and 5; M. Giavalisco et al., in preparation), while the size of our U dropout samples 
are relatively small due to the shallow depth of our ground-based data sets. We plan to ex- 
tend our observation of clustering segregation to the higher-luminosity regime at the same 
redshift intervals (z ~ 3 - 5) in a follow-up paper (K.-S. Lee et al., in preparation). 



-17- 



5.2. Small scale statistics 

Our deep ACS data also allow us to probe the existence of sub-structure in the halos 
that host the most massive LBGs. A general feature of hierarchical clustering is that more 
massive halos have higher probability than smaller ones to host more than one galaxy (e.g. 
Berlind & Weinberg 2002). If the existence of the clustering segregation is due to the fact 
that LBGs with higher UV luminosity are, on average, more massive than fainter ones, 
then one prediction is that brighter LBGs should also show a tendency to have more closely 
associated galaxies than fainter ones. Our deep ACS sample seems well suited to carry out 
this test, because satellite galaxies are typically less massive than the central galaxy of the 
halo (Berlind & Weinberg 2002), and because less massive galaxies apparently have fainter 
UV luminosity. It is likely that pairs, or multiplets, in which one or more members are 
significantly fainter than the brightest one are under-represented in the shallower ground- 
based samples. Unfortunately, the lack of systematic redshift measures prevents us from 
identifying individual physical pairs. Our measure will have to be statistical in character, 
and signal-to-noise will be lost due to projection. 

In our samples, pair statistics on small scales reflect not only large-scale structure of 
dark matter (via projected close-pairs) but also the physics of galaxy formation within given 
halos (via physical close-pairs). In this section, we discuss small scale statistics of the ACS 
selected samples only (-B435 and V606~band dropouts). 

As shown in Figure 2, a significant departure from a power-law on small scales is 
observed in our w(9) measures, most pronounced at the scale of 9 < 10 arcsec. 10 arcsec 
corresponds to roughly 0.3 - 0.4 hr x Mpc comoving at redshift 3-5. Recently, Ouchi et al. 
(2005b) detected a similar small-scale excess in their w(9) measurement of z ~ 4 galaxies 
selected from a deep wide-field survey in SXDF field (Ouchi et al. 2004). Steepening on small 
scales seems to be present in all samples, which we verify by fitting the w(9) measurement 
of each sample to a power-law when the data points with separations 9 < 10 arcsec are 
included, in which case we find that the slope is substantially steeper. The slopes of the 
power-law fits in both cases are shown in Figure 4 for the -B435 and Vgo6 _ band dropouts. 

This is consistent with the theoretical expectation that there are two separate contri- 
butions to w(9) or £(r), namely a one-halo term and a two-halo term, the origins of which 
are fundamentally different in nature. The two-halo term reflects the spatial distribution 
of host dark halos by counting galaxy pairs that belong to two distinct halos. This term 
dominates large-scale behavior of the correlation function and vanishes on small scales due 
to the halo exclusion effect, namely two halos cannot coexist closer than their typical size. In 
contrast, the one-halo term accounts for galaxy pairs sharing the same halo and is dominant 
on small scales. This reflects how galaxies populate halos. In this context, the single power- 
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law nature of the correlation function seen in the local universe is somewhat coincidental 
because it requires a smooth transition between these two terms to result in pure power-law 
which is possible only for a certain distribution of galaxies therein. Berlind & Weinberg 
(2002) argued that the contribution of galaxy pairs in massive galaxy clusters is critical for 
a smooth transition between these two terms. Clusters must be rare at high redshift, hence 
this steepening should be more prominent. This is also observed in high-resolution N-body 
simulations. For example, Kravtsov et al. (2004) measured the correlation function of ha- 
los above a given mass threshold identified from such simulations at various cosmic epochs. 
They found that at z = the halo correlation function can be fairly well approximated by 
a power law at all probed scales. However, at higher redshifts, steepening occurs on small 
scales, at progressively smaller scales with higher redshifts. As a result, power-law fits using 
the range of scales 0.1-8 h^ 1 Mpc give systematically steeper values of (3 compared with 
the fits over the range ~ 0.3 - 8ft- -1 Mpc. In fact, Figure 4 shows a similar behavior (in the 
angular domain) as Figure 12 of Kravtsov et al. (2004). 

Although typical halos at high-z must be much less massive than their present-day 
counterparts (therefore harboring fewer galaxies per halo), halos at high-mass end may host 
multiple galaxies. This is consistent with the recent discoveries that protoclusters were 
already forming by redshift of 4 - 6 (e.g. Ouchi et al. 2005a; Miley et al. 2004). To test for 
this possibility, we have looked for close galaxy pairs that may share host halos. In light of 
the luminosity segregation, we have split the full sample based on the Z850~band magnitude 
because statistically selecting bright galaxies is equivalent to selecting the sites of the most 
massive dark halos. Since the spectroscopic information is not available, we count galaxy 
neighbors around the brightest LBGs as a function of angular separation only. 

We split the galaxies in the full ACS sample (z ~ 4) into three sub-groups, z S50 < 24.3, 
m* < z S50 < 25.0 and m* + 0.5 < z S50 < 27.0 where m* for the -B 435 -band dropouts is roughly 
24.7 (M. Giavalisco et al, in preparation). The number of galaxies in each group is 29, 93 and 
2150, respectively. We consider galaxies in the two bright groups as candidates for central 
galaxies in halos, and the ones in the faintest group as satellite galaxies. We count the 
number of faint neighbors (z 850 > m * + 0-5 ~ 25.2) around the first two bright groups as a 
function of angular separation. The upper left panel of Figure 5 shows the average number 
of neighbors (cumulative counts) around these two bright groups up to 50 arcsec with la 
Poisson errors, however, these errors should be considered strictly as lower limits because 
these are correlated. What is expected for an uncorrected population with the same surface 
density is also shown with the dashed line. Due to the clustering, both z 850 < 24.3 and 
24.7 < z$5o< 25.0 are well above the dashed line, however, z$5o< 24.3 galaxies have more 
neighbors than the other (^ 2.5a). This cannot be understood in terms of the difference 
in their (large-scale) clustering strength because the same galaxies (therefore with the same 
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clustering strength) are counted as faint neighbors. In addition, there is clearly no reason 
to believe that bright galaxies have more galaxies along the line of sight which would be 
included as projected close-pairs. Another way to confirm that the difference does not come 
from clustering is to look at the same pair counts on larger scales. The upper right panel 
of Figure 5 shows the counts at 100 - 160 arcsec where the two slowly converge. This is 
because on large scales the contribution from halo-halo clustering dominates the counts thus 
reflecting the overall clustering strength. To remove the effect coming from the clustering, 
we take the ratio of the two which is shown in the lower panels. The ratio exceeds unity at 
9 < 30 arcsec and slowly goes back to unity on larger scales. 

This is compelling evidence that we are indeed observing central-satellite galaxy pairs 
sharing the most massive halos. This is also a natural outcome of the luminosity segregation 
in that the most massive halos (hosting the brightest LBGs) are likely to have sub-structures 
massive enough to form satellite galaxies bright enough for us to detect. We have tried a 
similar neighborhood search for the [/-band dropouts, but failed to detect the same trend. 
Because the ?7-band dropouts have larger area coverage thus larger number of bright sources, 
we have increased the magnitude cut for the brightest group up to R = 22.0 but still the 
brightest group does not evidently have more neighbors than the rest. It is possible that 
the relatively shallow depth of our [/-band data plays an important part in the lack of 
faint neighbors. If, for simplicity, no evolution of the luminosity function is assumed from 
z ~ 3 to z ~ 4, one needs [/-band dropouts fainter than R = to* + 0.5 = 25.0 (to* ~ 24.5; 
Adelberger & Steidel 2000) to carry out neighbor counts at a similar luminosity level to 
the i?43 5 -band dropouts case. In other words, we may be seeing only the bright end of the 
luminosity function where there is essentially one-to-one correspondence between host halos 
and galaxies. It is consistent with Adelberger et al. (2005) 's data that have comparable 
depth based on which they concluded that there is little evidence of galaxy multiplicity 
and that the fraction of dark halos hosting more than one galaxy is about 5%. Another 
possibility, though minor, is that the relatively poor seeing of our R band, which is where 
galaxy detection is carried out, prevents detection and deblending of pairs that are closer 
than 2 arcsec. 



6. The Halo Occupation Distribution at z ~ 4 and 5 

So far we have discussed the measures of the angular correlation function of LBGs, 
the implied spatial correlation length, and we have reported evidence of a departure of the 
correlation function from a single power law. We have also shown direct evidence of sub- 
structure in the most massive halos in the form of number counts of faint LBGs around 
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the brightest ones on scales 6 < 10 arcsec. Finally, we have also measured the clustering 
segregation, namely the fact that brighter LBGs are more strongly clustered in space. In 
order to understand these results in light of underlying physics, we need to compare with 
theoretical predictions through which quantitative interpretations of the result can be made. 

Associating the spatial distribution of galaxies with that of dark halos requires an under- 
standing of the relationship between galaxy properties and those of host halos, for example, 
galaxy colors, stellar mass or luminosity are correlated with circular velocity, total mass of 
halos and so forth. We do not yet know if all the halos above a certain mass threshold host 
galaxies with similar SEDs (colors) and luminosity to those of LBGs to be included in our 
sample. At z ~ 3, a population of red galaxies, either heavily obscured by dust or old, 
dubbed as DRGs, which are bright in the near-infrared have been found at similar redshifts 
with about half the space density of [/-band dropouts (1Z < 25.5) (Franx et al. 2003; van 
Dokkum et al. 2003). Many of these galaxies would not have been selected via the Lyman 
Break technique. Moreover, the submillimeter galaxies, though their space density is much 
smaller ((2.4 ± 1.2) x 10~ 6 Mpc~ 3 ; Chapman et al. 2003), also coexist in the same redshift 
range. These high-z populations (including LBGs) suffer from respective selection biases 
thus cannot be solely representative of high-z population in general. The halo occupation 
distribution (HOD) is a simple model that sets a statistical relationship between halo mass 
and its probability to form galaxies, thus provides us with a powerful tool to understand the 
characteristics of different populations independently and study their respective relations to 
host halos. The HOD formalism has been successfully applied to various local galaxy samples 
(e.g. different luminosities, colors, spectral types etc.; Zehavi et al. 2004b, Magliocchetti et 
al. 2003) and also at intermediate redshift (QSOs; Porciani et al. 2004). Note that the HOD 
formalism can be extended to the conditional luminosity function (CLF; Yang et al. 2003, 
van den Bosch et al. 2003), to understand the relationship between halo mass and galaxy 
luminosity. We are currently working to constrain this relationship at z ~ 3,4 and 5 using 
very large samples of LBGs extracted from the COSMOS 2-square degree survey, which 
provide excellent sampling at the high end of the luminosity (and hence mass) distribution 
(K.-S. Lee et al. in preparation). 

We move on now to consider a simple model for the halo occupation distribution func- 
tion, which can be used to understand the phenomenology of the observed clustering prop- 
erties of LBGs in physical terms, and derive its parameters. These include the minimum 
mass of halos that can host the observed galaxies and the average halo mass for each galaxy 
sub-sample. We will then use the best-fit HOD parameters to discuss the inferred halo bias 
and its implications. 
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6.1. The HOD Formalism and its Physical Implications 



The angular correlation function can be considered to have two separate contributions, 
galaxy pairs from the same halo and from two distinct halos. 



On large scales, only W2h{d) contributes to the total angular correlation function because any 
galaxy pair with a large angular separation cannot reside in the same halo. On the other 
hand, on small scales and on scales comparable to virial diameter of halos, both terms will 
be present with a ratio that depends on the degree of galaxy multiplicity. The two-halo term 
can be constrained from clustering measurements; however, the galaxy multiplicity function 
or halo occupation distribution (HOD) is not very well constrained at high redshifts because 
of the still relatively small and shallow sample. Bullock et al. (2002) estimated the HOD 
from a sample of 800 LBGs with spectroscopic redshifts at z ~ 3 based on a ground-based 
survey (7Z < 25.5; Steidel et al. 1998). Although this sample has the advantage of having 
spectroscopic information, the median luminosity of the sample is significantly brighter than 
our ACS sample at z ~ 4, and thus the number of faint members of pairs or multiplets is 
very likely quite small. Adelberger et al. (2005) estimate that at most 5% of their LBGs 
at z ~ 3 (an extended version of the sample used by Bullock et al. 2002) are members of 
close pairs. We do not yet know the luminosity function of satellite LBGs; however, our 
neighborhood analysis discussed earlier for the -E?435-band dropout sample at z ~ 4 showed 
that a large fraction of the excess number of neighbors comes from -2g5o> 26 (galaxies counted 
as faint neighbors have 25.2 <z 850 < 27.0). Using our template LBG spectrum, ^850 = 26.0 
corresponds to TZ = 25.5 at z ~ 3, that is the limiting flux of Adelberger et al.'s sample. 
While the bright members of pairs and multiplets are certainly present in their sample, 
Adelberger et al.'s sample must be highly incomplete for satellite galaxies. 

Once the halo power spectrum is given and an HOD model assigned, one can unam- 
biguously determine the correlation function (Seljak 2000; Berlind & Weinberg 2002; Berlind 
et al. 2003) with suitable assumptions as to how galaxies are distributed within halos (see 
later). We adopt a simple functional form (Jing et al. 1998; Wechsler et al. 2001). This 
model assumes that all galaxies are associated with dark matter halos and makes a simpli- 
fying assumption that the number of galaxies within a given halo depends only on the halo 
mass and can be modelled by a single power-law: 



w(6)=w lh (9)+w 2h (6) 



(6) 



(M/Mi) a if M>M, 



mm 



= 



otherwise 



(7) 



where Mi is the typical halo mass at which there is on average one galaxy per halo, M min is 
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the cutoff halo mass below which halos cannot host galaxies and a is the power-law index, 
the efficiency of galaxy formation within halos of given mass M. 

While these three parameters, M 1: M min and a are sufficient to constrain the correlation 
function on large scales, the shape of the correlation function on small scales also depends on 
how galaxies populate host dark halos. This is because the one-halo term w\h is proportional 
to (Ng(N g — 1)) while W2h is proportional to (N g ) (see Appendix A. for details). Berlind 
& Weinberg (2002) discussed what these effects are and their impact on the shape of the 
correlation function on small scales. This includes the scatter of N g around the mean, 
galaxy concentration with respect to halo concentration, the existence of a central galaxy at 
all times and the galaxy velocity dispersion with respect to that of halos. These factors can 
be used to construct more sophisticated models to constrain the HOD once a larger galaxy 
sample becomes available. However, at high redshift, the number of galaxies included in any 
sample (so far) remains much smaller than that in such surveys as 2dF or SDSS, hence a 
simple approach seems preferable. Therefore we make simplifying assumptions to this model 
as follows. We force the first galaxy to be placed at the center of host halos and the rest 
(satellite galaxies) largely follow the halo mass distribution for which we adopt an NFW 
profile (Navarro, Frenk & White 1997). The latter is consistent with the finding of Kravtsov 
et al. (2004) that the sub-halos in their simulations follow the distribution of dark matter 
within the halos. 

In order to calculate w(6) we closely follow the recipe given by Hamana et al. (2004), 
using the procedure described in Appendix A. We generate a library of model w(8)'s for a 
wide range of (Mi, a) and determine the best-fit w{6) that matches the abundance (n g ) and 
the clustering strength (M min , a and Mi) simultaneously using x 2 minimization. Note that 
M min is not a free parameter once the space density is specified. 



where N(z) is the redshift selection function obtained from the simulations and n h (M) is 
the halo mass function for which we adopt the analytic halo model by Sheth & Tormen 
(1999, 2002). Note that we derive the HOD parameters, M x and a, from the w(9) of the full 
sample (B270 and V270) and do not attempt to constrain these parameters separately for 
the sub-samples. We also fix the observed number density to the incompleteness-corrected 
number density integrated down to z 850 = 27.0 for the fit. 

Table 4 lists the HOD parameters obtained using this method for the B 435 and V^-band 
dropouts only. For the [/-band dropouts, we are unable to put any robust constraints on 





/ dzN(z)[dV(z)/dz]n g (z) 
JdzN{z)[dV{z)/dz] 



(8) 
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these parameters primarily because the HOD fit is most sensitive to w(9) at small separations. 
The w{6) measure of the ?7-band dropouts on such scales seems to indicate a constant value 
rather than steepening. This is more consistent with the correlation function having the 
two-halo term only, i.e. each observed galaxy belongs to a different halo, in which case 
W2h stays relatively constant on small scales actually falls off but the projection, w(6), 
stays constant due to the contribution from projected galaxy pairs). Figure 6 shows the best- 
fit model w{6) for the full sample of B i35 and V6 6~band dropouts with the observational 
measures. The corresponding confidence levels for the HOD fits are illustrated in Figure 7. 

The best-fit parameters for the full samples are (M m j„, a, Mi) = (7x 10 w h~ 1 M Q , 0.65, 
l.SxfCP/i^M©) for the 5 43 5-band dropouts and (5 x KP/i^Mq, 0.80, f.O x lO^/i^M©) 
for the V606~band dropouts, respectively. Though galaxies were selected with different filter 
systems and also their sample is shallower, it is interesting to compare them to those of 
Hamana et al. (2004), (1.6 x lO 11 /^ 1 !!©, 0.5, 2.4 x IO^/i^Mq), (1.4 x Kpfc^M©, 0.5, 
1.4 x 10 12 /i~ 1 M Q ) for their LBGz4s and LBGz5s samples. M min is found to be a factor 
of 2 - 3 smaller from our survey. This is hardly surprising because M m i n depends mainly 
on the depth of given survey. Our sample is significantly deeper than theirs (~ 1 mag) 
and, because LBGs have compact morphology, more complete in the magnitude range where 
there is overlap. As a result, the overall clustering strength of our sample is weaker which 
is most directly reflected in M min . The parameter Mi, the mass scale at which statistically 
there is one galaxy corresponding to each halo, is smaller in our case, consistent with our 
expectations. On the other hand, this work and previous studies consistently show that a 
should be 0.5 -1 though it may be dependent on M (thus the depth). Therefore over a 
reasonable range of observable halo mass it can be approximated (to the zeroth order) as 
a constant, thus nearly independent of mass scale. Unfortunately, as can be seen in Figure 
7, it is the least constrained parameter from the HOD fit. In principle, however, a can be 
constrained independently via studies such as our neighborhood analysis or gravitational 
lensing surveys. 

We calculate (N 9 )m (note the subscript "M" to distinguish the quantity from the pre- 
viously mentioned (N g (M))), the ratio of galaxy number density to halo number density 
i=n g /nh) integrated over the allowed range of halo mass. We also calculate (M), the aver- 
age mass of halos that host the observed LBGs (assuming the best- fit HOD parameters) as 
follows: 



$M mi M M )^ M )) dM 



g/M — poo 



fM mi M M ) dM 

J-_M(N 9 (M))n h (M)dM 
{ } f~(N g (M))n h (M)dM W 
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For each sub-sample, we compute M min that matches the number density, then {N 9 )m 
and (M) are computed using the best-fit HOD parameters. Table 4 summarizes the results. 
Uncertainties are computed by propagating the bootstrap errors of w(6) after marginalizing 
over Mi and a. We find (N 9 )m = 0.3 and 0.2 for the full sample of 5 4 3 5 -band and Vqqq- 
band dropouts, respectively. If each galaxy belongs to a different halo, this quantity would 
represent the average number of galaxies per halo. For example, for our S 435 -band dropouts, 
this means that only 30% of the halos harbor LBGs bright enough to be detected (in other 
words, 30% of the halos that are more massive than M min ). In reality, where multiple galaxies 
can share host halos, {N g ) M can be regarded as an upper limit. Note that Adelberger et al. 
(2005) and Giavalisco & Dickinson (2001), based on the ground-based sample with 1Z < 25.5, 
argued that most halos at z ~ 3 host visible LBGs. While it is not clear if their claim is in 
quantitative discrepancy with our estimate, we note that possible evolutionary effects from 
z ~ 4 to z ~ 3 remain poorly understood at this time. It is interesting, however, that the 
i?435-band dropouts with z^q< 26.0 (B260), which have a comparable absolute luminosity 
threshold to that of the [/-band dropouts of 1Z < 25.5 (U255), have (N 9 )m ~ 0.5, similar 
to what Adelberger et al. (2005) have found (see their Figure 10). This suggests that LBGs 
at z ~ 3 and 4 of similar luminosity also have a similar duration of star-formation during 
which they are visible (duty-cycle hereafter). Note that if the average galaxy luminosity is 
closely tied to the halo mass, as suggested by the clustering segregation, then in principle one 
can empirically find the intrinsic relationship between UV luminosity and mass at a given 
redshift 2 . For example, given a mass spectrum and a HOD function, such a relationship 
will simulataneously reproduce the shape of the correlation function, that of the clustering 
segregation, and the luminosity function. Interestingly, when the same absolute magnitude 
cut is made to the V<306~band dropouts, we find that (N 9 )m is a factor of 2 - 3 smaller. 
We postpone the interpretations of these results to the next subsection. Finally, we note 
explicitly that the clustering properties of LBGs are consistent with the notion, built in the 
HOD model itself, that (N g ) M becomes smaller for less massive halos. A possible physical 
explanation for this fact could be that in low-mass halos, gas cooling (hence star formation) 
is much more inefficient due to photo-ionization (squelching) or SNe feedback (e.g. Silk 
1997, Efstathiou 2000). In addition, it is possible that the duty-cycle of galaxies may be 
mass-dependent such that lower-mass halos have shorter duty-cycle, thus resulting in a 
smaller value of (N 9 )m- 



2 This would be the average "observed" UV luminosity, not corrected for the effects of dust obscuration. 
Since not all galaxies have the same dust obscuration, and since the amount of dust obscuration is, in general, 
dependent on the UV luminosity itself, the statistical relationship between obscured and unobscured UV 
luminosity and mass will be, in general, different. 
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6.2. A Different Approach: Separation of Central and Satellite Galaxies 

We have also explored an alternative approach, which is to model central galaxies and 
satellites separately (Guzik & Seljak 2002; Kravtsov et al. 2004; Zheng et al. 2005) In this 
model, the probability for a halo to have a central galaxy is a step function, namely either 
the halo has a central galaxy or it does not, whereas the distribution of satellite galaxies 
follows a power-law: 

(N g (M)) = l + (M/M[f ifM>M min 

= otherwise (10) 

We refer to this model as the "one-plus" model as opposed to the previous one which we 
refer to as the "pure power-law" model hereafter. Note that the normalization M[ and 
the threshold mass M' min have different physical meanings in this model compared to the 
previous one. In the present model the parameter M[ is the halo mass at which there is 
one satellite galaxy (thus two galaxies in total) in the host halo. The parameter M' min is 
conceptually the same as M m i n of the first model, except that now halos above this threshold 
value always have (N g ) > 1, while in the pure power-law model only a fraction of halos with 
M > Mi are host to more than one galaxy on average. Using the alternative model for the 
halo distribution function of Lyman-break galaxies we find that the best-fit parameters for 
the 5435-band dropouts are {M' min , a 1 , M[) = (2 x 1O 11 ^ 1 M , 1.3,4 x lO^/i^M©). 

We have fit the two HOD models to our data, and the \ 2 analysis shows that both 
reproduce the observed w(6) equally well. Figure 8 shows the predicted w{6) from the 
best-fit HOD parameters in both cases together with the observed correlation function of 
the -B 4 3 5 -band dropouts. We found xl+ — 0-97 an d X 2 vv i — 0-91. In addition, they predict 
similar (N g ) for mass scales of 10 12 - 10 13 h~ 1 M Q , to which most of massive halos included 
in the sample belong. This is shown in the right panel of Figure 8. At the high-mass 
end (M > 1O 13 /i -1 M ), however, the two models diverge from each other rather quickly. 
To constrain the HOD function in this mass range one needs a much larger sample covering 
significantly more area than the GOODS, because of the rarity of high-mass halos. On scales 
M < 10 12 h~ 1 M G , on the other hand, the "pure power-law" model predicts numerous low- 
mass halos being included in the sample with a relatively low efficiency of actually hosting 
a galaxy, while for the "one-plus" model, {N g ) is truncated (by definition) at the mass scale 
of M < 2 x I0 11 h- 1 M Q . 

Finally, we observe that although the two models are equally good at reproducing the 
observations, this is very likely the result of the still relatively large uncertainties in the 
measure of w(9). The rationale behind the "one-plus" model is that every single halo above 
M min participates in forming observed galaxies. In contrast, the "pure power-law" model 
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allows a fraction of halos to form galaxies over some range of mass (M < Mi), with the 
complementary fraction remaining quiscent (either because no star formation has not started 
in them or becuase observed in between bursts of star formation). The currently available 
data at high rdshifts do not allow us to test which model is a better representation of reality. 
However, we note that the built-in assumption of the "one-plus" model, that there is always 
one central galaxy (in our case, one LBG), in every halo above a mass threshold, naturally 
precludes the possibility of having other types of galaxies (i.e. non-LBG) as the central 
galaxy in those halos. While such a limitation is not a problem in the local universe, it is 
very likely a misrepresentation of the situation in the high-redshift universe. For example, 
a number of observations have shown the existence of massive galaxies at z ~ 3 whose UV 
spectral energy distribution is such that they would not be selected as Lyman-break galaxies. 
Such galaxies very likely have masses that are in the high-mass end of the LBG mass range 
(e.g. Shapley et al. 2004, van Dokkum et al. 2004). The "one-plus" model would not 
provide a fair representation of the HOD if such galaxies are relatively common at high 
redshift. For these reasons we will limit the discussion only to the "pure power-law" model 
in what follows. 

Realistically, we expect observed LBGs to be the central galaxies in only a certain 
fraction of halos (M > Mi). For this reason, we expect that (N g ) = p+(M/Mi) a (M > M min 
with < p < 1) may be a more reasonable description of the actual halo occupation 
distribution such that a fraction 1 — p is reserved for other populations. While M 1 and 
a (to a lesser extent) may depend on the depth of a given survey, p more likely represents 
a fundamental quantity that is the efficiency of LBGs as central galaxies on the high end 
of the halo mass function. We plan to explore this possibility in a follow-up paper with a 
larger data set (K.-S. Lee et al., in preparation). 



6.3. The Evolution of Halo Bias with Redshift 

We compute the nonlinear autocorrelation function of matter £ m (r) at redshift z ~ 3, 
4 and 5 by adopting the algorithm by Peacock & Dodds (1996). £ m (r) is then inverted to 
w m (6) using the Limber transform as follows: 

, a , rdzN\z)JZ [dx/R H (z)]U[DM(z)er + xY /2 
Wm{9) = [jdzN{z)Y (U) 



where Dm(z) is the proper motion distance and Rh(z) the Hubble radius at given redshift 
z. Using this method, the average bias is calculated as a function of angular separations. 
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w(6) 



(12) 



w m {9) 



The bias (b b s (9)) monotonically decreases at small separations then stays constant at 
9 > 70 arcsec. We adopt the value of (b obs ) at 100 arcsec as our linear bias value. The 
theoretical predictions of the galaxy and halo bias can be computed from the given halo 
model as follows: 



where n h = J™ mm n h (M)dM and n g = J^ m , n n h (M)(N g (M))dM. We use the best-fit HOD 
parameters for (N g (M)). Note that halo bias (bh) is HOD-independent whereas (b g ) is not. 
The observed galaxy number density is calculated by fitting the full sample to the Schechter 
luminosity function after correcting for photometric incompleteness (e.g. Giavalisco 2005) 
and integrating down to each magnitude cut (see Table 3). The error bars in the number 
density are derived by propagating the uncertainties from the Schechter LF fit. 

Figure 9 shows the average galaxy bias from the clustering measures together with the 
predictions of (b h ) and (b g ) at z ~ 4 and 5. Note that the abscissa represents the halo 
abundance (rih) for the dashed lines and the galaxy abundance (n g ) for the solid lines and 
symbols (n g — nh if N g = 1). The observed galaxy bias values are in good agreement with 
the CDM model predictions when the best-fit HOD parameters are assumed. 

To compare the results from different cosmic epochs, it is important to remove the 
effect of the luminosity segregation which could be mistakenly interpreted as "evolution". 
We truncate each sample to the same absolute luminosity threshold to match that of our 
brightest sample, the [/-band dropouts (R < 25.5), the absolute magnitude of which is 
^1700 < -20.0. This corresponds to R < 25.5, z 850 < 26.0 and z 850 < 26.6 at z ~ 3, 4 and 5, 
respectively. Since these luminosity thresholds are similar to those of the U255, B260 and 
V265 sample, we compare the bias values for these magnitude cuts. Figure 10 illustrates how 
these bias values change with redshift compared to halo biases for various mass thresholds. 
The dashed lines indicate the evolution of the average halo bias with a constant halo mass 
threshold for M halo = 10 10 ,5 x 10 10 , 10 n ,5 x 10 11 , 10 12 and 5 x 10 12 h~ x M Q (from bottom) 
from the Sheth & Tormen (1999) model. Lines serve as as a guide to determine the effective 
halo mass threshold of a given galaxy sample. 



(b g ) = — / n h (M)b h (M){N g (M))dM 



1 f°° 



g Jm, 




(13) 



It is intriguing that our measures at z ~ 3 and 4 imply that these galaxies are hosted by 
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halos of mass 5 x lO n /i _1 M < M ha i Q < lO 12 /i _1 M , however, at z ~ 5, the host halo mass is 
approximately lO n /i _1 M , a factor of 5 - 10 times smaller than its lower- z counterparts at a 
given fixed luminosity. Hamana et al. (2004) reported the average bias of their LBGz4s with 
similar absolute luminosity threshold, i' < 26.0, to be 3 - 4.5, consistent with our results. 
This implies that at z ~ 5, star-formation may have been more efficient, hence galaxies 
of comparable luminosities were hosted in much lower-mass halos. This result is, in fact, 
another manifestation of what we have mentioned in an earlier subsection, that (N g ) was a 
factor of 2 - 3 smaller at higher redshift when the same absolute luminosity cut was made 
for a given choice of HOD parameters. 



7. Summary 

1. We have studied the spatial clustering properties of LBGs at three different cosmic 
epochs. We used three samples of LBGs which include two deep ACS-selected LBGs samples 
at z ~ 4 and 5 down to a magnitude limit of -285o~ 27 and the ground-based sample at z ~ 3 
down to R ~ 25.5. The ACS samples enable us to study the spatial distribution of the 
large-scale structure as well as of small-scale clustering whereas the ground-based sample 
allows us to probe the large-scale clustering of bright LBGs selected from a large contiguous 
area. 

2. We found that the clustering strength (the amplitude of w(9) or value of r ) of LBGs 
at all redshifts that we have studied depends on the UV luminosity (rest-frame A ~ 1700 A) 
of the galaxies. This luminosity segregation is in good quantitative agreement with previous 
investigations based on much brighter galaxy samples and we are able to robustly extend 
this reesult down to at least 1 mag fainter than previous studies. The faintest galaxy sub- 
sample has a real-space correlation length as small as 2 h~ l Mpc, only half of that known 
for 1Z < 25.5 sample. A physical implication is that galaxies brighter at UV wavelengths are 
hosted by more massive dark halos, implying that the star-formation is primarily regulated 
by local gravity, and suggesting that other physical mechanisms (major mergers, interactions) 
are secondary drivers. 

3. We have carried out neighbor counts around bright LBGs in our sample and counted 
the number of bright-faint galaxy pairs. We discovered that the likelihood of finding faint 
neighbor galaxies (z^q> m* + 0.5) around bright ones (z^q< 24.3) up to scales comparable to 
the virial radius of dark halos, is at least 20 - 30 % higher than random LBG-LBG pairs with 
a high significance. This is consistent with the interpretation that there are faint compan- 
ions around the brightest galaxies sharing host halos, supporting that the departure of the 
correlation function on small scales (from the extrapolation of the large-scale power-law) is 
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indeed due to central-satellite pairs. The observed sub-structure is also consistent with the 
luminosity segregation that we have detected, in that brighter galaxies (more massive halos) 
are likely to have fainter companions (halo sub-structures). We do not yet know the rela- 
tionship between the UV luminosity and the actual host halo mass. However, the clustering 
segregation gives support to the interpretation that the most luminous LBGs are hosted by 
the most massive halos and tend to have more satellite galaxies in halo sub-structures. Due 
to the relatively small dynamic range in magnitude of our survey, we are limited to a bright 
cut of -2850= 26.0 and also to cumulative magnitude thresholds for our samples, however, we 
expect this effect to be more dramatic in future surveys that cover larger cosmic volume, 
which should confirm our results. 

4. Our measurement of the correlation function at z ~ 4 and 5 indicates a steepening 
on small scales, thus cannot be well described by a simple power-law. We also find that the 
scale at which the steepening occurs is comparable to the angular size of dark halos at the 
relevant redshift range, and hence we attribute this steepening to the presence of multiple 
galaxies within the same halo. A simple HOD model with a suitable scaling law between halo 
mass and the number of galaxies seems to be able to explain both the large-scale and small- 
scale behavior of the angular correlation function. We find the best-fit HOD parameters for 
LBGs at z ~ 4 and 5 to be (M min , a, Mi) = (7 x l0 10 h~ 1 M Q , 0.65, 1.3 x lO^/i^M©) and 
(5 x lO lo /i _1 M , 0.8, 1.0 x 1O 12 /i _1 M ), respectively. With these parameters, we find that, 
on average, the number density of galaxies is a factor of 4 - 5 smaller than that of dark halos 
for the full sample and that this discrepancy is smaller for brighter galaxy samples. This 
suggests that the duty-cycle of faint LBG populations may be shorter than that of brighter 
counterparts suggested by shallower surveys, but at the same absolute luminosity cut, our 
results are in good agreement with shallower surveys. However, at z ~ 3, our results do not 
indicate a similar steepening, but is more consistent with a constant value for the correlation 
function at small scales. We interpret this discrepancy to be caused by the difference in 
survey depth, as shallower surveys tend to pick up only bright LBGs which probably tend to 
be central galaxies and not satellites (which are usually fainter) within halo sub-structures 
that largely produce this small-scale feature. 

5. The scaling of the correlation length or galaxy bias as a function of volume density 
shows a similar trend to that expected for dark halos, implying that our measures fully 
support the biased galaxy formation scenario predicted by the cold dark matter framework. 
This suggests that LBGs flag the sites of dark halos efficiently, however, not every halo is lit 
up at the same time. This has to do with the fact that the duty-cycle of LBGs is shorter 
than the cosmic time spanned by the observations (0.3 - 0.5 Gyr), therefore only a fraction 
of halos are lit up by star-formation at any given time. In addition, a HOD that scales 
with halo mass might be an indication that duty-cycle is also mass-dependent such that 
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higher-mass halos may have longer duty-cycles than their lower-mass counterparts. 

6. We truncate each sample to a fixed absolute magnitude and compute the average 
bias of these galaxies. At z ~ 5, galaxies with a given luminosity are hosted by halos of mass 
~ lO u h~ 1 M Q whereas at z ~ 3 and z ~ 4, the corresponding host halos are a factor of 5 - 
10 more massive. The implication is that at z ~ 5 star-formation was more efficient than 
the later epochs. 

We thank Cristiano Porciani for very helpful assistance. KL acknowledges Takashi 
Hamana for kindly providing his halo code for model comparisons. MG and KL thank 
David Weinberg for his very useful comments on HOD modelling. We would also like to 
thank an anynymous referee for useful suggestions. We are grateful to the entire GOODS 
team and particularly, Tomas Dahlen, who has put great efforts into calibrating many of the 
GOODS ground-based data sets, which made our ground-based measures possible. 



A. w(9) and Halo Occupation Distribution 

For a given set of (Mi, a) with a fixed number density (n g ), the minimum mass M min is 
determined by matching the number density of the observed galaxies (n g ) with that expected 
from the halo model as follows: 

POD 

n g {z) = / dMn h (M)(N g (M)) (Al) 

_ J dzN(z)[dV/dz]n 9 (z) 

JdzN(z)[dV/dz] [ ] 

where N(z) is the normalized redshift selection function and dV/dz is the comoving volume 
element per unit solid angle. The contributions from the one-halo and two-halo components 
are computed from their respective galaxy power spectra, P g h (k,z) and P g h (k,z). w{6) is 
the inverse Fourier transform of the total power spectrum, P g (k,z): 

w(6) = J dzN(z) 2 ^)- 1 J dk^-P g (k,z)J (r(z)6k) (A3) 

where Jo is the Bessel function of the first kind and r(z) is the radial comoving distance. 
The one-halo contribution to the power spectrum is given as: 

Pl\Kz) = — I dM^(N g (N g -l))\y(k,M)\v (A4) 
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where y is the normalized Fourier transform of the halo mass profile (Seljak 2000; Scherrer 
& Bertschinger 1991). We further assume that each halo has an NFW profile with the 
concentration parameter varying with the halo mass as found by Bullock et al. (2001). The 
parameter p is taken to be p = 1 if (N g (N g — 1)) < 1, and p = 2 otherwise (Seljak 2000). 

The two-halo term contribution to the galaxy power spectrum according to the linear 
halo bias model (Cole & Kaiser 1989; Mo & White 1996) is: 

P* h (k, z) = Pun(k){^j J dM^(N 9 (M))b(M)y(k, M)} 2 (A5) 

where Pu n (k) is the linear dark matter power spectrum and b(M) is the linear halo bias from 
the fitting function of Sheth & Tormen (1999). 
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Redshift 

Fig. I. — The redshift distribution functions N(z) of the U, -B435 and V^-band dropouts 
(from top) estimated from Monte Carlo simulations. The coarser histograms indicate the 
distribution of spectroscopically confirmed objects to date whose numbers are specified on top 
right corner of each panel. The histograms from the simulations are arbitrarily normalized 
to match the distribution of the spectroscopic samples. 
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Fig. 2. — The angular correlation function w(0) of the B i35 (the left panel) and V^-band 
dropouts (the right panel). The points are the measurements corrected for the integral 
constraint (IC) together with the best-fit power-laws. The solid line shows the results when 
the slope is allowed to vary and the dashed line shows the result when the slope is fixed 
to 0.6. The fit is done including angular separations 6 > 10 arcsec only. While the data 
can be well described by a single power-law on large scales, a significant departure from a 
power-law on scales 9 < 10 arcsec is observed. 
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Fig. 3. — Left panels: The angular correlation function w{9) of the U, -B435, V606~band 
dropouts (from top). The filled symbols show the measures from the full sample (ACS- 
based: zg5o< 27.0, ground-based: R < 25.5) while the open symbols show results from the 
brightest sub-sample (ACS: -2 850 < 26.0, ground-based: R < 24.0) together with the best-fit 
power-law when the slope is forced to be (3 = 0.6. The correlation amplitude, A w , increases 
by a factor of more than 2 for the brighter groups in each case. Right panels: The real-space 
correlation length r is shown as a function of the magnitude thresholds used to define each 
sample. Note that the top-right panel is i?-based (R < 24.0, 24.5 and 25.5), whereas the 
others are z 8 5o-based (z 8 5o< 26.0, 26.5 and 27.0). 
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Fig. 4. — The steepening of the correlation slope on small scales: both panels show the best- 
fit slope of the correlation function for samples with a given apparent magnitude threshold. 
The filled (open) symbols show the slope (3 when fitted to a single power-law without (with) 
the observational measures at angular separations of 9 < 10 arcsec. This is due to an 
excess number of physically close galaxy-pairs causing a significant departure from a pure 
power-law at small separations. 
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Fig. 5. — The top panels show the average number of faint LBG neighbors (-2850> m* + 0.5; 
m* = 24.7 at z ~ 4) as a function of angular separation around the two bright groups, 
z 850 < 24.3 (filled) and m* < z 850 < 25.0 (open), together with what is expected for purely 
random distributions (w — 0), shown as the dashed lines. The error bars represent Poisson 
errors. The bottom panels show the ratio of the two quantities shown in the top panels. 
Taking the ratio of the two essentially removes the effect from the projected galaxy pairs, 
providing a better indication of the true excess due to physical pairs. The galaxies that 
belong to the brightest group (-2850< 24.3), on average, have more faint LBG neighbors 
(3c significance) than the other group (m* < z$5o< 25.0). This effect is most pronounced 
at < 6 < 30 arcsec. The ratio slowly drops down to unity on larger scales. This is 
consistent with the halo sub-structure interpretation because when the angular separation 
is comparable to the typical distance between nearby halos, (N neighbor ) is dominated by 
large-scale clustering rather than halo sub-structure. 
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Fig. 6. — The two-component fit with the best-fit HOD parameters is shown together with 
the observed angular correlation function. The dashed line shows the one-halo contribution 
and the dashed-dot line shows the two-halo contribution to the final w(9). The solid line is 
the sum of the two contributions. The best-fit HOD parameters are M min = 7 x lO lo /i _1 M , 
Mi = 1.3 x lO^/i^M© and a = 0.65. Right: the V 60e -band dropouts: The best-fit HOD 
parameters are M min = 5 x lO^/i^M©, M x = 1.0 x lQ 12 h' l M & and a = 0.80. 
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Fig. 7. — The x 2 contour map of the HOD fits overlaid with confidence intervals for the -B435 
(left) and V(}06 _ band (right) dropouts: Mi vs. a. The contour lines show A% 2 = 2.3 and 6.2, 
corresponding to 68% and 95% confidence levels. Boxes on the bottom right corners show 
the best-fit HOD parameter values and the galaxy number density. 
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Fig. 8. — Left panel: The two-component fit of the "pure power-law" (thick lines) and "one- 
plus" (thin lines) model is shown together with the observed angular correlation function of 
the 5435-band dropouts. The two models describe the observed correlation function equally 
well (x?+ — 0.97 and Xp P i — 0.91). Right panel: the mean occupation number (N g ) predicted 
from the two best-fit HOD models is shown as a function of mass. Halo mass is in units of 
h~ x M Q . The two models predict similar (N g ) values in the mass range of 10 12 - 10 13 h~ l M Q) 
to which most of the massive halos "observed" in the data most likely belong. This implies 
that our estimate is very robust regardless of the employed HOD model in the high-end 
of the halo mass function. However, at the mass scale of M > iO 13 /i _1 M , (N g ) from the 
"one-plus" model largely exceeds the other. On the other hand, at M < 10 11 /i" 1 M o , the 
"pure power-law" model implies a larger number of low-mass halos in the samples, while 
the "one-plus" model is truncated (by definition; see text for further discussion). 
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Fig. 9. — The average bias vs the galaxy (or halo) abundance: The data points and the solid 
lines represent the average galaxy bias from the correlation measures and the predictions as a 
function of galaxy number density, respectively. The average bias is computed as linear bias 
(Sheth & Tormen 1999) weighted by the best-fit HOD to account for the galaxy multiplicity. 
The dashed lines show the average halo b function of halo number density. If there 

were one-to-one correspondence between halos and galaxies, then the two would coincide. 
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Fig. 10. — The evolution of the average bias as a function of redshift: the data points 
indicate the average bias implied for the U, -B435 and V606 _ band dropouts at the fixed absolute 
luminosity threshold M 1700 < —20.0, which corresponds to 1Z < 25.5 at z ~ 3. The dashed 
lines show the predicted linear bias from the Sheth & Tormen (1999) model for halo mass 
M ha i > M min = 10 10 ,5 x 10 10 , 10 n ,5 x 10 11 , 10 12 and 5 x 10 12 K^Mq (from bottom). At 
z ~ 3 and 4, the observed galaxies are hosted by dark halos of similar mass (5 x 10 11 h~ 1 M Q 
- 10 12 /i _1 M ). At z ~ 5, the observed clustering strength of the galaxies implies that they 
are hosted by less massive halos, M haio ~ IO^^Mq. This suggests that star-formation 
may have been more efficient at z ~ 5. 
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Table 1: The space-based data set 



Field 




A a 


^(5 4 35) b 




^77 5 ) b 


a-(z 850 ) h N LBG / 


NlBG,5 C 


GOODS- 


-N 


160.0 


29.00 


29.08 


28.33 


28.09 1169 


502 


GOODS 


-S 


160.0 


28.88 


29.05 


28.34 


28.06 1294 


376 



a Survey area, in units of arcmin 2 

b l a surface brightness fluctuations within 1 arcsec 2 apertures 
c Number of i?43 5 /V606 _ b an( i dropout Lyman-break galaxies 



Table 2: The ground-based data set 



Field A a(U) a a(B) a cr{R) a N v h 

extended GOODS-S 1200 29.42 29.75 28.95 1609 



a l a surface brightness fluctuations within 1 arcsec 2 apertures 
h Number of [/-band dropout Lyman-break galaxies (R < 25.5) 
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Table 3: The LBG abundance, the angular correlation function and the inferred correlation 
lengths 



Flavor Sample n 9 (xlO 3 ) A W)0 (3 r 00 A W;1 (3\ r 0i i 



u 


R < 


25. 


.5 


3. 


,3± 


1. 





0. 


p; fi +0.04 
JU -0.04 





.6 


4. 


1 +0.1 

■-■--0.2 


0. 


PT9+0.28 
JZ -0.20 


0. 


6 o+0.11 

.uo_ 12 


4 


q+0.2 


u 


R < 


24. 


,5 


0, 


,7± 


0. 


.3 


1 


1 a+om 

iu -0.08 





.6 


6 


c+0.2 
°-0.3 


0. 


r/i+0.23 
°^-0.16 





44+0.O8 

•^-0.08 


6 


4+0.3 

■^-0.3 


u 


R < 


24. 





0. 


,2± 





.1 


1 


cro+0.14 
•°°-0.17 





.6 


7. 


O+0.4 
•O-0.5 


0. 


o/f+0.68 
° 4 -0.32 


0. 


c.9+0.12 
oz -0.11 


7. 


O+0.5 
• 8 -0.6 


-B435 


^850< 


27. 


.0 


7. 


,3± 


1 


.0 





qo+0.04 
OO -0.05 





.6 


2 


q+0.2 

y -0.2 


0. 


49+0.32 
^ z -0.26 


0. 


fiO+0.16 
uy -0.15 


2 


o+0.2 
°-0.2 


-B435 


^850< 


26. 


,5 


4. 


,5± 





,7 





fif) +0.10 
uu -0.06 





.6 


3 


O+0.3 
y -0.3 


0. 


7^+0.56 
' J -0.37 





7n +0.17 
' ' u -0.13 


3 


7+0.3 
• ' -0.3 


-B435 


^850< 


26. 


.0 


2. 


,7± 





.4 





qq+o.14 

yy -0.13 


0. 


.6 


5. 


o+0.4 
•°-0.5 


0. 


64+0-98 
°^-0.45 


0. 


•°^-0.19 


5. 


1+0.4 
• x — 0.5 


^606 


^850 < 


27. 


.0 


4. 


,2± 





.8 





• ' °-0.15 





.6 


4. 


4 +0.5 
•^-0.5 


1. 


nn+1.68 
uu -0.75 





oc+0.20 
■ oo -0.23 


4. 


9+0.4 
• z -0.5 


^606 


^850 < 


26. 


,5 


2. 


6± 


0. 


.5 


1. 


1 9+0.34 
• iz -0.25 





.6 


5 


o+l.l 
•O-0.8 


2. 


n o+8.80 
■ U °-1.76 


1 


1 0+030 
• iu -0.27 


4. 


4+0.7 

• 4 -0.7 


^606 


^850 < 


26. 





1 


,5± 


0. 


.3 


1 


7n +0.42 
• ' u -0.37 





.6 


7. 


■°-1.0 


3.* 


jo+12.89 
^-3.46 


1 


1 0+0.31 


5 


o+l.l 
■°-1.0 



Note. — The number density (n g ) and the comoving correlation lengths (r ) arc in units of h 3 Mpc~ 3 
and h" 1 Mpc, respectively. The quantities with the secondary subscript "0" are when the slope is fixed to 
a fiducial value (3 = 0.6, those with the subscript "1" are when both (3 and A w are allowed to vary. 



Table 4: The best-fit HOD parameters a , the ratio of galaxy number density to halo number 
density {N g ) M , and the average halo mass (M) b for different luminosity thresholds 



Flavor 


Sample 


M ± h ' C 


a c 


(N g ) M 


(M) b 


-B435 


^850 < 27.0 


1.3xl0 12 


0.65 


0.31±0.14 


(4±l)xlO n 


-B435 


^850< 26.5 


1.3xl0 12 


0.65 


0.38±0.13 


(6±l)xlO n 


-B435 


^850< 26.0 


1.3xl0 12 


0.65 


0.49±0.13 


(8±2)xlO n 


^606 


z 850 < 27.0 


l.OxlO 12 


0.80 


0.20±0.18 


(2±l)xlO n 


^606 


^850< 26.5 


l.OxlO 12 


0.80 


0.24±0.19 


(3±l)xlO n 


^606 


^850< 26.0 


l.OxlO 12 


0.80 


0.30±0.21 


(4±2)xlO n 



"Using the "pure power-law" model 
b All masses are in units of /i~ 1 M Q 

c Note that the same HOD parameters, Mi and a, constrained from the full sample, were used for all sub- 
samples 



